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In this article we show how Ehrenfest mean field theory can be made both a more accurate 
and efficient method to treat nonadiabatic quantum dynamics by combining it with the generalized 
quantum master equation framework. The resulting mean field generalized quantum master equation 
(MF-GQME) approach is a non-perturbative and non-Markovian theory to treat open quantum 
systems without any restrictions on the form of the Hamiltonian that it can be applied to. By 
studying relaxation dynamics in a wide range of dynamical regimes, typical of charge and energy 
transfer, we show that MF-GQME provides a much higher accuracy than a direct application of 
mean field theory. In addition, these increases in accuracy are accompanied by computational speed- 
ups of between one and two orders of magnitude that become larger as the system becomes more 
nonadiabatic. This combination of quantum-classical theory and master equation techniques thus 
makes it possible to obtain the accuracy of much more computationally expensive approaches at a 
cost lower than even mean field dynamics, providing the ability to treat the quantum dynamics of 
atomistic condensed phase systems for long times. 

PACS numbers: 


I. INTRODUCTION 

The exact treatment of real time nonadiabatic quan¬ 
tum dynamics in condensed phase chemical systems re¬ 
mains a significant challenge that spurs the ongoing de¬ 
velopment of approximate methods that are accurate, ef¬ 
ficient, and can treat systems with a wide range of differ¬ 
ent forms of interactions. In particular quantum-classical 
(semiclassical) trajectory based methods offer a hierarchy 
of approaches, derived from the exact real time path inte¬ 
gral formulation of quantum mechanics, that offer differ¬ 
ent balances between accuracy and computational cost. 

At the lowest tier of this hierarchy is Ehrenfest mean 
field theory (MFT), which neglects all dynamical corre¬ 
lations between between the quantum (subsystem) and 
classical (bath) degrees of freedom [T| . Above this lie lin¬ 
earized path integral approaches such as LSC-IVR. im 
FK-LPI [J], PBME [5j, and LAND-map JBj which, al¬ 
though still mean field in nature, capture some correla¬ 
tion. These methods offer higher accuracy than mean 
field theory, at the expense of at least an order of mag¬ 
nitude more computational effort arising from the need 
to average over the mapping variables associated with 
the quantum subsystem degrees of freedom. Recently it 
has been shown that partially linearizing the propagator 
for the electronic degrees of freedom, giving rise to the 
partially linearized density matrix (PLDM) [7 approach 
and the forward-backward trajectory solution (FBTS) to 
the quantum-classical Liouville equation (8j, allows more 
dynamical correlation to be included at relatively little 
additional cost to fully linearized methods. Introduc- 
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ing further dynamical correlations between the subsys¬ 
tem and the bath adds further accuracy, at the expense 
of assigning weights and phase factors to the trajectories 
HM2I- This in turn adds many orders of magnitude to 
the number of trajectories, which grows rapidly with sys¬ 
tem dimensionality and time, that must be generated in 
order to obtain converged properties. 

Since moving up this hierarchy requires orders of mag¬ 
nitude more computational effort, only the lowest tiers 
are likely to be practical, both now and in the foresee¬ 
able future, for nonadiabatic problems containing large 
quantum subsystems or where on-tlie-fly treatment of the 
electronic states is required. At present, this means that 
one is typically limited to using MFT, or Tully’s fewest 
switches surface hopping (FSSH) algorithm and its vari¬ 
ants pEHIZ], or linearized path integral approaches fl8l - 

m- 

The generalized quantum master equation (GQME) of¬ 
fers an alternative way of exactly describing the nona¬ 
diabatic evolution of a quantum subsystem by formally 
recasting the effects of the environment into a memory 
kernel [211 [22]. In the condensed phase the environment 
comprises of a large number of modes with a broad range 
of frequencies that couple to the quantum subsystem, 
leading to a memory kernel that typically decays much 
more rapidly than the population relaxation time of the 
subsystem. This separation in time-scales becomes more 
pronounced as the system-bath coupling strength or the 
nonadiabaticity is increased. In such regimes, trajectory 
based approaches suffer from the rapid accumulation er¬ 
rors in lower tier methods and a rapid rise in the com¬ 
putational cost with time in higher tiered approaches. 
Hence, using trajectory based approaches to calculate 
the memory kernel, which is short lived compared to the 
subsystem relaxation time, and then using the GQME to 
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generate the subsystem dynamics, offers massive advan¬ 
tages in terms of both accuracy and computational cost 
compared to a direct application of the trajectory based 
approach. Such a realization has previously been shown 
to be highly effective in the case of higher-tier trajectory 
based methods l23ff2o1 . However, it also allows one to 
make the most of the lower tier approaches by allowing 
a more accurate treatment of strongly nonadiabatic and 
coupled problems while retaining their existing strengths 
in treating weakly coupled and adiabatic regimes. 

Here we show that MFT can be combined with the 
GQME formalism to yield a method that is as accu¬ 
rate as higher tiered trajectory-based techniques such as 
FBTS and PLDM, but requires a much lower computa¬ 
tional effort than MFT alone, without being limited to 
any particular form of the Hamiltonian. We show how 
to calculate the memory kernel of the GQME using dy¬ 
namical trajectories obtained by solving the mean field 
equations of motion for the system, and how to use these 
mean field kernels to subsequently generate the reduced 
dynamics of the quantum subsystem. Finally we demon¬ 
strate that our method is more computationally efficient 
and more physically accurate for treating charge and en¬ 
ergy transfer regimes of the spin-boson model than a di¬ 
rect application of MFT. 


II. THEORY 

A. Ehrenfest Mean Field Theory 

A particularly simple and instructive route to derive 
the Ehrenfest mean field equations of motion is to begin 
with the quantum-classical Liouville equation |26j and 
neglect correlations in the system-bath dynamics. One 
begins by considering a system in which only a small 
subset of the degrees of freedom behave quantum me¬ 
chanically, and are of interest. This set of degrees of 
freedom is denoted as the quantum subsystem (or sim¬ 
ply as the subsystem), and the remainder of the system 
is referred to as the bath, and is assumed to behave es¬ 
sentially classically. 

The total Hamiltonian for the entire system is written 
as a sum of subsystem, bath, and coupling terms, 

H = H s + H b + H sb , (1) 


The quantum-classical Liouville equation [26| , 

-pw(X,t) = -iCp w (X,t), (3) 

describes the time evolution of the density matrix 
pw ( X , t ), which is a quantum mechanical operator that 
depends on the classical phase space variables. The 
quantum-classical Liouville (QCL) operator is 

'] — "} — {') H w }), (4) 

where [•,•] is the commutator, and {•,•} is the Poisson 
bracket in the phase space of the environmental variables. 
The subscript W refers to the partial Wigner transform 
over the environmental degrees of freedom in the system. 
The partial Wigner transform of the density operator, p, 
is 

Pw(R,P) = / dZe iPZ (R-^\p\R+^)- (5) 

In order to arrive at MFT, one makes the approxima¬ 
tion that the density of the system can be written as a 
product of the subsystem and bath reduced densities at 
all times. This asserts that there are no dynamical corre¬ 
lations between the subsystem and the bath, as the total 
density remains in a product state. 


Pw(X, f) = p s (t)p b (X,t), (6) 


where the bath RDM is pb(X,t) = Tr a (pw{X,t)). Using 
the approximations in Eqs. (3) and (6), one obtains the 
Ehrenfest mean-field equations of motion m 


dp s {t) 
dt 



dXH sb (R(t))p b (X,t) \ ,p s (t) , 


dR{t) _ P(t) 
dt M ’ 


dP = d{H b , w {X(t)) + Tr s {H sb {R{t))p s {t)}) 
dt dR[t ) 


The mean field approximation for a subsystem observ¬ 
able, O s (f), can be written as 


(O s (t)) = J dX (Tr s {O s p s (t)}) Pb{X,t). (8) 


where the subscripts s, b , and sb refer to the subsystem, 
the bath, and the system-bath coupling, respectively. 
The time-evolution of the reduced density matrix of the 
subsystem is defined as 

p s (t) =Tr b (p(t)) = J dXp w (X,t), (2) 

where p is the density operator for the entire system 
and pw is its Wigner transform, Tr b indicates the par¬ 
tial trace taken over the bath degrees of freedom, X = 
( R,p ) = (-Rl, i?2) PN b , Pi, P 2 , PN b ): and N b is the 
number of bath degrees of freedom. 


MFT is a highly efficient approximation to quantum dy¬ 
namics which, while obtaining accurate results in some 
regimes, fails when quantum effects in the bath are im¬ 
portant and in regimes with a nonzero subsystem energy 
bias. 

B. The Generalized Quantum Master Equation 

An alternative approach to formulating the expecta¬ 
tion value of a subsystem observable in terms of the den¬ 
sity matrix for the full system, is to treat the reduced 
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dynamics of the subsystem via projection operator tech¬ 
niques. One begins with the exact quantum evolution of 
the full system, governed by the Liouville - von Neumann 
equation, 

= (9) 

and focuses only on the evolution of the subsystem de¬ 
grees of freedom by projecting out the degrees of freedom 
of the bath. This can be accomplished by applying a pro¬ 
jection operator, V , 

V = p e h q ®Tr b {-) = pl« w f dX{-), (10) 

and its compliment, Q = 1 — V. 

The form of the coupling part of the Hamiltonian, H sb 
is chosen to be 

H sb = S® A, (11) 

where S' is a pure subsystem operator, and A is a pure 

bath operator. We use this factorization for simplicity 
and, since any coupling Hamiltonian can be written as a 
sum of such operators, there is no loss of generality [28] . 
We also write the bath part of the coupling Hamiltonian, 
A, such that its equilibrium thermal average vanishes, 

(A) eq =Tr b [Ap e b q ] = 0, (12) 

In the problems that we consider in Sec. m such a con¬ 
dition is naturally satisfied, but in general this condition 
can be enforced by redefining A relative to its thermal 
average [29] . 

We restrict our considerations to initial states of the 
Feynman-Vernon type, where the initial density of full 
system can be factorized in the following manner, 

P(t = 0) = p s (0) ®p e b q , (13) 

where 

= exp (~pH b ) 

Tr b [exp(-/3H b )\ 

is the density operator for the isolated bath in thermal 
equilibrium. Under these conditions the exact time evo¬ 
lution of the subsystem RDM is given by the Nakajima- 
Zwanzig GQME [23 [22] 

d A 

—p s {t) = —iC s p s (t) — / dTK{r)p s {t - t), (15) 

at J o 

where the subsystem Liouville operator, C s , is given by 
C s = \[H S , •], and the memory kernel, JC, is given by 

K,(t) = Tr b {C sb exp {-iQCT)QC sb p e b q } (16) 

The assumption of an initially uncorrelated state is not 
necessary, and relaxing it would introduce an inhomoge¬ 
neous term to the GQME which depends on the chosen 


initial state (that can also be calculated by mean field 
theory and other approximate methods). The free sub¬ 
system evolution prescribed by C s is generally simple to 
simulate, and hence calculating the evolution of the sub¬ 
system RDM reduces to the calculation of the memory 
kernel, which encodes the deviation from free subsystem 
evolution. From Eq. (15), it is clear that changes in the 
subsystem populations in the GQME are driven by the 
subsystem Liouville operator as well as the memory ker¬ 
nel and thus the lifetime of the population dynamics is 
typically longer than that of the memory kernel. In con¬ 
densed phase systems, where the environment spans a 
wide spectral bandwidth, the memory kernel is expected 
to decay much more quickly than the population relax¬ 
ation time. 

The general form for the memory kernel, given above, 
is not straightforward to evaluate since it explicitly de¬ 
pends on the projection operator. However, there are a 
variety of ways that it can be constructed from simula¬ 
tions of the unprojected dynamics of the full system. We 
will limit our discussion of this procedure to a particular 
method introduced by Shi and Geva, who showed that 
the full memory kernel, /C(t), can be written in terms of 
a set of partial memory kernels [23 : . 1281 ; 59], 

IC(t ) = /Ci(t)-H f dr , /Ci(r - r')/C 2 (r), (17) 

Jo 

/C 2 (r) = /C 3 (t)-H [ dr'/C 3 (r - r')/C 2 (r), (18) 

Jo 

where the partial memory kernels are given by 

Mr) = Tr b {C sb e~ iCT C sb p eq }, (19) 

IC 3 (r) = Tr b {e~ iCr C sb p e b q }. (20) 

By combining Eq. (15), and Eqs. (17 - 20) the subsys¬ 
tem RDM can thus be exactly evolved using projection- 
free input. Exact numerical evaluations of these expres¬ 
sions have recently been carried out using techniques such 
as QUAPI, real-time quantum Monte Carlo, and ML- 
MCTDH [59H35]. 

C. Mean Field Evaluation of the Memory Kernel 

In practice, evolving the subsystem RDM using the 
GQME formalism for an arbitrary system is no less cum¬ 
bersome than solving Eq. (9). One way to proceed, that 
is valid in the weak coupling limit, is to factorize the 
memory kernel into subsystem and bath parts which can 
be evaluated separately, leading to the well established 
Bloch-Redfield theory [25] [33] . In this limit IC 2 and K ,3 
vanish, and the memory kernel can be evaluated using 
equilibrium correlation functions of the isolated bath. In 
contrast, here we make no such simplifying assumption 
and instead evaluate the kernels using MFT. 

Since K, 1 and /C 3 do not contain any projected input, 
they can be simulated directly and then used to obtain 
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/C by solving Eqs. (17) and (18). The matrix elements 
of /Ci and /C 3 are obtained by projecting each quantity 
onto a basis which spans the subsystem Hilbert space, 

(/CiUw'(r) = (5^(T)A£'(r)5^(0)A(0)} e9 
-(^ Q ,(r)A^'(r)^(0)A(0)) e9 
+ (K(0)S^(0)k^{T)S^a'(T)) eq 
-(A(0)S^ M (0)A^(r)5 a ^(r)) e? , 

( 21 ) 

(JC 3 )aa'^(r) = ((l b )^“'(r)^(0)A(0)} e? 

-(5^(0)A(0)(i 6 )^'(r)) e 9 ,(22) 

where a, a', (3, and j3' refer to subsystem states, and If, 
is the unit operator for the bath. In the above two ex¬ 
pressions the Einstein summation convention is used. 

The matrix elements of the partial memory kernels, 
/Ci and /C 3 , contain correlation functions of the following 
form, 


(S(0)A(0)fJ'(r)) e , = Tr^H sb \p)(p'\ 

xe tCT/h f\a')(a\y (23) 

where T is a bath operator (If, or A). Defining operators, 
A = H s b(h ® \P)(P'\) and B = i& < 8 > \a'){a\ , or B = 
H s b(ib < 8 > |a')(a|), expression (23) takes on the general 
form for a quantum time correl 


lation function, 
(5(0)A(0)f5'(r)) eg = Tr{ffAB{T)). 


(24) 


where the equilibrium density corresponds to that of the 
isolated bath. 

Working in the coordinate representation of the bath 
degrees of freedom, and making use of the partial Wigner 
transform, Eq. (24) can be rewritten as 


2 . K -2 is generated from /C3 by an iterative solution 
to Eq. (18), using /C3 itself as an initial guess for 
/C 2 . This iterative procedure typically converges 
very quickly, and often requires only a few tens of 
iterations. 


3 . /Ci and /C2 are used as input to obtain the full mem¬ 
ory kernel /C by numerical integration of Eqs. (17) 
and (18). 


4. Using the full memory kernel, the evolution of the 
subsystem density is generated by direct numerical 
integration of the GQME using Eq. (15). 


In our calculations the kernel elements were calculated 
for the specified time r, and set to zero for all t > r; no 
smoothing of the kernel data was performed for t < r. 
Using the MF-GQME approach one can propagate the 
subsystem RDM for arbitrarily long times using only 
short-time information obtained from mean field trajec¬ 
tories. 


III. RESULTS AND DISCUSSION 

In order to assess the accuracy and efficiency of our 
MF-GQME approach, we performed simulations of the 
spin-boson model. Despite its apparent simplicity, this 
system is a prototypical model for the study of quan¬ 
tum transport and relaxation processes in the condensed 
phase O 125], and remains a challenging test to ap¬ 
proximate methods. Since it is now possible to gener¬ 
ate numerically exact results in many of the parameter 
regimes of the spin-boson model, it provides an ideal 
benchmark test case for the accuracy and efficiency of 
approximate nonadiabatic dynamics approaches. In par¬ 
ticular we compare our MF-GQME approach to a direct 
MFT treatment, as well as to the recently introduced 
FBTS method j8j [9, which has been shown to outperform 
fully linearized methods at marginal extra computational 
cost. 

The spin-boson Hamiltonian can be written in the sub¬ 
system basis as 


Tr(pl 9 AB(t)) = Tr s J dX [p e b q A] w (X, 0 )B W (X, r)(25) 

The calculation of the memory kernel then amounts 
to the evaluation of the above expression, which can be 
performed by a hybrid Monte Carlo / molecular dynam¬ 
ics algorithm where (i) initial conditions are sampled 
from \pl Q A] w {X, 0), and (ii) the system is propagated 
in time from these initial conditions using MFT to evalu¬ 
ate Bw(X,t). The full memory kernel can then be con¬ 
structed, and the subsystem RDM propagated as follows: 

1. Mean-field trajectories are used to obtain the corre¬ 
lation functions necessary to form /C 1 and /C 3 using 
Eqs. (21) and (22). See Appendix A for more de¬ 
tails on the calculation of these quantities for the 
model studied here. 


H — e<j z + A cr x + A? c jRj&z \ 1 

j ' ' 

( 26 ) 

where a x and a z are Pauli spin matrices. This Hamil¬ 
tonian describes a two level quantum system with ener¬ 
getic bias 2 e, and electronic coupling (tunneling) matrix 
element A, that is bi-linearly coupled to a bath of inde¬ 
pendent harmonic oscillators. In the spin-boson model 
both subsystem states are coupled to the same bath in 
an anti-correlated fashion. In contrast to problems where 
the subsystem states are coupled to independent uncor¬ 
related baths, such as the Frenkel exciton model, this 
form of the coupling presents a more difficult challenge 
for mean-field type methods to describe lilMEB- The 
greater challenge of treating anti-correlated baths is due 
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Figure 1: Evolution of the subsystem population difference 
in the biased, nonadiabatic regime with increasing system- 
bath coupling strength; lj c = 2A, t = A, /3 = 5A _1 . In 
each panel the exact results are shown in solid black dots, the 
dotted blue lines are direct MFT results, the solid red lines 
are FBTS results, and the solid blue lines are the results of 
our MF-GQME approach. 


to the greater difference between the mean force and 
those on the individual diabatic surfaces compared to 
independent, uncorrelated, baths. 

The interaction between the system and the bath can 
be fully characterized by the spectral density, J(w), 
which determines the strength of the interactions be¬ 
tween the subsystem and bath, which we chose to be 
of Ohmic form, 


JM = (27) 

The Kondo parameter, £, controls the strength of the 
coupling between subsystem and the bath, and the cut¬ 
off frequency co c sets the primary time-scale for the bath 
evolution. The quantum subsystem was initialized in di¬ 
abatic state 1 , and the bath was sampled from its (iso¬ 
lated) equilibrium distribution. In our calculations 400 
bath modes were used to represent the continuous spec¬ 
tral density which, for all regimes and approaches em¬ 
ployed, gave results converged to graphical accuracy. 

Figure [T] compares the subsystem evolution as the 
system-bath coupling is increased, for a system with 
an energetic bias (e 7 ^ 0 ) in the nonadiabatic regime 
(^ > 1). Numerically exact results were generated us¬ 
ing our own implementation of the quasi-adiabatic path- 
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Figure 2: Matrix elements of the memory kernel of the GQME 
for for u>c = 2A, e = A, /3 = 5A _1 , £ = 0.4, 5 = 0.02A -1 , and 
Ntraj = 2xl0 4 . In each panel the MFT results are shown as 
blue lines, and the exact QUAPI results are the black lines. 
Note the different y-axis scales in each panel, due to the vary¬ 
ing magnitudes of each element. 


integral (QUAPI) algorithm [42H44] , The performance 
of trajectory based approaches degrades as the system 
becomes more nonadiabatic and the system-bath cou¬ 
pling is increased, a failure that is particularly evident 
for systems with nonzero energetic bias. In this regime 
direct MFT is over-coherent and only captures the exact 
QUAPI results at very short times (r A < 1.5). The long 
time population difference is close to zero, which is a no¬ 
torious feature of MFT that arises due to the deviation of 
the mean force from the forces on the individual surfaces, 
resulting in an accumulation of error in the subsystem 
population distribution as time progresses. FBTS, which 
is also mean field in nature but includes more dynamical 
correlation between the subsystem and bath, much more 
accurately captures the long time populations at a mod¬ 
est increase in the number of trajectories of by a factor 
of « 5 in this regime. In contrast to MFT, the FBTS 
underestimates the oscillatory nature of the population 
decay. 

By using MFT to approximate the memory kernel, 
our MF-GQME approach produces populations dynam¬ 
ics that are in perfect quantitative agreement with the 
exact quantum results. Figure [2] shows the MFT mem¬ 
ory kernel elements that give rise to the dynamics in the 
bottom panel of Fig. [TJ compared with the exact QUAPI 
results. There are four nonzero, linearly independent el¬ 
ements of K, for a two-level system (due to the form of 
the spin-boson Hamiltonian, K, aa ppi =0). As expected, 
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Figure 3: Evolution of the subsystem population difference in 
the intermediate coupling regime with increasing bath nona- 
diabaticity; e = A, /3 = 5A _1 , £ = 0.1. In each panel the 
exact results are the solid black dots, the dotted blue lines 
are direct MFT results, the solid red lines are FBTS results, 
and the solid blue lines are the results of our MF-GQME ap¬ 
proach. 


MFT fails to correctly capture the long time behavior; 
exhibiting spurious oscillations when tA > 2. This is 
consistent with the fact that the direct MFT treatment 
of the population dynamics in Fig. [I] also begins to show 
marked errors at those times. Despite the deviations of 
MFT from the exact kernels, the population dynamics 
generated using a memory kernel of length tA = 1.5 
(shown in the bottom panel of Fig. [T]) are in excellent 
agreement with the exact results. Using memory times 
longer than tA = 1.5 includes the spurious long time os¬ 
cillations but only introduces very minor changes to the 
population dynamics. As MF-GQME only requires the 
generation of very short trajectories to obtain the mem¬ 
ory kernel, the entire population decay, which occurs in 
a time of approximately 15A -1 , can be obtained at a 
cost 10 times cheaper than a standard mean-field calcu¬ 
lation of the same observable and ~ 50 times cheaper 
than FBTS. 

The effect of increasing the characteristic frequency of 
the bath, w c , which pushes the system into an increas¬ 
ingly nonadiabatic regime, is displayed in Fig. [3j In the 
adiabatic limit MFT is accurate, but as the nonadiabatic- 
ity is increased the direct MFT results begin to deviate 
from the exact solution at progressively shorter times, 
and by w c = 7.5A it is unable to even capture the first 
minimum in the coherent decay correctly. FBTS again 
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Figure 4: Evolution of the subsystem population difference 
generated from the MF-GQME approach using memory ker¬ 
nels of varying length; e = A, /3 = 5A _1 , £ = 0.1. 


performs better than direct MFT in reproducing the long 
time limit, albeit at an increased cost of an order of mag¬ 
nitude more trajectories, but again gives results which 
are increasingly overdamped as w c is increased. This 
overdamping is consistent with that observed in previ¬ 
ous FBTS and PLDM results in similar regimes pirn. 

The MF-GQME results are in quantitative agreement 
at low nonadiabaticity, and even at the highest nonadi- 
abaticity exhibit only a very subtle phase shift relative 
to the exact results. Again this reflects that the mem¬ 
ory kernel decays rapidly in these nonadiabatic regimes, 
and hence retaining a memory kernel of length r A = 1.5 
is sufficient to generate the results shown in Fig. [3] 
Figure [4] shows the convergence of the population dy¬ 
namics obtained from the MF-GQME approach when 
different amounts of time, r A, are included in the mem¬ 
ory kernel for the regimes shown in the top and bottom 
panels of Fig. [3] In both cases the convergence is essen¬ 
tially monotonic as the length of memory in the kernel 
is increased. The MF-GQME results are generally bet¬ 
ter than direct MFT even for very short memory kernels. 
When an insufficient length of time is used in the ker¬ 
nel, the error accumulated in the propagation of the sub¬ 
system RDM in MF-GQME manifests in the observed 
population dynamics more prominently at longer times. 

For the more adiabatic regime shown (top panel) the 
dissipation induced by the bath is well captured when 
r A = 0.7, whereas when the system becomes more nona¬ 
diabatic (bottom panel) this is decreased further with 
convergence obtained around rA = 0.5, which is 30 
times shorter than the population decay time. This 
again highlights the complementarity of using MFT and 
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Figure 5: Evolution of the subsystem population difference 
at zero temperature in the unbiased case, e = 0. In each 
panel the exact ML-MCTDH results from Ref |45j are the 
solid black dots, the dotted blue lines are direct MFT results, 
the solid red lines are FBTS results, and the solid blue lines 
are the results of our MF-GQME approach. 
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Figure 6: Marcus electron transfer regime, with uj c = 10A, 
£ = 1.0, /3A = 0.05, and A = 2£u} c is the reorganization 
energy. The top panel shows the equilibrium subsystem pop¬ 
ulation difference after relaxation has occurred while the bot¬ 
tom panel shows the electron transfer rate as a function of 
driving force obtained from exponential fits to the population 
decay. In the top panel the solid black line is the Boltzmann 
distribution, in the bottom panel the solid black line is the 
Marcus rate, the dotted blue lines are MFT results, the solid 
red circles are FBTS results, and the solid blue squares are 
the results of our MF-GQME approach. 


other trajectory based approaches in conjunction with 
the GQME framework; moving further into the nonadia- 
batic regime, which leads to a faster breakdown of MFT, 
requires less total time to be included in the GQME mem¬ 
ory kernel. 

In contrast to the spin-boson model at finite tempera¬ 
ture, the zero-temperature limit is considered to be more 
challenging, as nuclear quantum effects in the bath be¬ 
come more prominent; a classical treatment of the bath 
distribution would predict only a single allowed bath ini¬ 
tial configuration. In all the results shown here, the 
Wigner sampling of the bath ensures that zero-point en¬ 
ergy is included exactly in the initial condition, although 
this is not guaranteed to be preserved by the approx¬ 
imate evolution prescribed by MFT or FBTS. In Fig. 
[5] we compare to exact multi-layer multi-configurational 
time-dependent Hartree (ML-MCTDH) results at zero- 
temperature in the nonadiabatic regime [35] ■ As seen in 
Fig. 0 and Fig. [ 3 ] increasing the nonadibaticity, <u c /A, 
or system-bath coupling, £, degrades the performance of 
MFT and FBTS at progressively shorter times although 
due to the lack of energetic bias the long time limits are 
in, probably fortuitously, good agreement. Suprisingly, 


FBTS performs worse than MFT in these regimes and 
again exhibits over-damped behavior resulting in slower 
relaxation. The MF-GQME results are in excellent agree¬ 
ment, with some very mild underdamping present in the 
case where the system-bath coupling is strong (middle 
panel). 

In addition to the significant increase in accuracy af¬ 
forded by MF-GQME over FBTS and direct MFT, this is 
accompanied by a significant increase in efficiency. As in 
Fig- 0 the memory kernel decay times in Fig. [5]required 
to obtain the converged population dynamics are sub¬ 
stantially shorter than the population decay times they 
generate, with memory kernel of lengths tA = 0.2, 0.2 
and 0.1 required to obtain the results in the top, mid¬ 
dle and bottom panels respectively. This again reflects 
that the memory kernel becomes shorter as the system 
becomes more non-adiabatic. Due to the much faster de¬ 
cay of the memory kernel than the populations dynamics 
(which occurs in approximately 5A _1 ), the MF-GQME 
results in Fig. 4 are between 25 and 50 times cheaper to 
generate than direct MFT dynamics and up to 500 times 
cheaper than FBTS. 

One of the most well-known difficulties of mean field 
theory occurs in the Marcus electron transfer regime 











146] . While most approximate dynamics approaches, like 
FSSH and MFT, are capable of qualitatively capturing 
the famous rate turnover as a function of driving force 
[1611251 071 05], as shown in the bottom panel of Fig. [6j 
MFT fails to correctly describe the donor-acceptor prod¬ 
uct ratios for the process. Much like the PBME approach 
[25], the FBTS method qualitatively captures the rate 
turnover but has a slightly asymmetric shape. Also, per¬ 
haps surprisingly given their mean field treatment of the 
system-bath interactions, both PBME and FBTS give 
the correct equilibrium distribution at long times. The 
ability of FBTS to qualitatively and semiquantitatively 
capture Marcus turnover is consistent with recent PLDM 
results for a similar model of Marcus electron transfer 
that was simulated using a flux-side correlation function 
approach [49] . The agreement between FBTS and PLDM 
is expected as the methods are identical when the sub¬ 
system Hamiltonian is chosen to be traceless [9]. They 
thus offer a similar tradeoff between cost and accuracy 
in the semiclassical hierarchy. 

While FBTS performs well by including some dynam¬ 
ical correlations between the system and bath, the com¬ 
plete neglect of these correlations in direct MFT leads to 
poor performance. This is strongly pronounced in biased 
regimes (those with large electronic driving forces) where 
the reaction rates are underestimated and the long time 
population difference is too small. The MF-GQME ap¬ 
proach resolves these issues, giving rise to quantitative 
agreement with the Marcus prediction for the rates and 
the Boltzmann distribution of the long time populations. 
In order to obtain the MF-GQME results in this regime 
r A = 0.5 time units of data were typically included in 
the memory kernels, and the population decay occurs on 
a timescale of a few hundred time units. This allowed for 
a ~ 200 fold efficiency gain compared to direct MFT, in 
addition to a substantially improved accuracy. 

IV. CONCLUSIONS 

Here we have shown that utilizing MFT within the 
GQME framework gives rise to nonadiabatic relaxation 
dynamics that are highly accurate across a wide range 
of physical regimes. This is accompanied by a computa¬ 
tional savings of one or two orders of magnitude com¬ 
pared with direct MFT calculations, and up to three 
orders of magnitude over FBTS. This success can be 
rationalized based on the fact that the subsystem Li- 
ouville operator is treated exactly within the GQME 
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VI. APPENDIX 


A. Explicit expressions for computing K ,i and K ,3 
for the spin-boson model 


According to the procedure outlined in Sec. (IIC), the 


integration in Eq. (25) is carried out over the phase space 
of the bath by Monte Carlo sampling initial conditions 
from [p h q A] w (X, 0) and generating dynamical trajecto¬ 
ries to evaluate Bw{X,t). 

If A is independent of momentum and linear in the 
coordinates of the bath, the Wigner transform of the op¬ 
erator product is [50] 


[Af>l q ]w = \.Pl q Aw = A w pl« 


dAw 


dpt q 


b,W 


2 H dR dP 


(28) 


which can be evaluated analytically for the harmonic 
bath employed in this study. 

In the spin-boson model, A = — JY CjRj and it’s 
Wigner transform is A^ = — JY Cj Rj. The Wigner 
transform of the equilibrium density for the isolated bath 
is 


= n 

j 


tanh (3u}j/2 

7 T 


x exp 


2 tanh /3ojj/2 


P 2 
2Mj 



(29) 


In practice, our initial conditions for the bath degrees 
of freedom were sampled from the initial bath density (by 
taking out a factor of the bath density from Eq. (28)) 
and the trajectories were then re-weighted using the re¬ 
maining term in Eq. (28). 

The subsystem operator in the system-bath coupling 
part of the Hamiltonian is S = a z , and the matrix ele¬ 
ments of S in the subsystem basis are given by 

(m\S\ri) = {m\a z \n) = S n S mn = (-l)" + 1 <5 mn . (30) 


where, m and n are eigenstates of the isolated subsystem, 
i.e. <r z |l) = | 1 ) and <r z | 2 ) = -| 2 ). 

The elements of K, 1 and K .3 are then given by 
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(ici) aa 'w(T) = ( { -ir« + (-ir^)fdx 

+{(-lf +a+1 + (-lf +a ') J dX 


n eq - 

i 

dA\y 

dPb!\V 

°b,W 

2 h 

Oil 

dP 

) eq 4 - 

i 

dAw 

®Pb,W 

6 ,W ^ 

2h 

OR 

OP 


and 


(K-'S'jaa' fi/3' ('7”) — 


dX 

A wpl q w - 

i dA w 

®Pbyv 

2 h dR 

dP 

dX 

A w P^w 

i dAw 

dpl q w 

2 h dR 

dP 


(X,0)A w (X,r)pf'(0)p“'«(r) 

(X,0)A w (X,r)pf'(0)p“' Q (T), 

(31) 

(X,0)pf'(0)p“'“(T) 
(X,0)pf'(0)p“'“(r), (32) 


where p““ = c a c* a , is the subsystem RDM, and c a and c a / are the coefficients of the subsystem wavefunction. 



